Atomic order of rare earth ions in a complex oxide: a path to magnetotaxial anisotropy

Complex oxides offer rich magnetic and electronic behavior intimately tied to the composition and arrangement of cations within the structure. Rare earth iron garnet films exhibit an anisotropy along the growth direction which has long been theorized to originate from the ordering of different cations on the same crystallographic site. Here, we directly demonstrate the three-dimensional ordering of rare earth ions in pulsed laser deposited (EuxTm1-x)3Fe5O12 garnet thin films using both atomically-resolved elemental mapping to visualize cation ordering and X-ray diffraction to detect the resulting order superlattice reflection. We quantify the resulting ordering-induced ‘magnetotaxial’ anisotropy as a function of Eu:Tm ratio using transport measurements, showing an overwhelmingly dominant contribution from magnetotaxial anisotropy that reaches 30 kJ m−3 for garnets with x = 0.5. Control of cation ordering on inequivalent sites provides a strategy to control matter on the atomic level and to engineer the magnetic properties of complex oxides.


Supplementary Information
This description of the site preference model has been summarized from the work of Herbert [3][4] (1) GIA occurs in complex oxides in which spin-orbit coupling (SOC) makes the magnetic properties sensitive to the type and arrangement of electrons in neighboring ions.Strong SOC also explains the existence of magnetocrystalline anisotropy. 5) GIA originates from the rare-earth ion sites (c sites) due to the different orientations of the coordination dodecahedra within the unit cell with respect to neighboring magnetic cations (Fe 3+ and RE 3+ ).At the growth surface, differently oriented c sites present different coordination constellations to incoming cations.Each unit cell has 24 c sites comprising 12 different orientations.The 12 types of sites can be obtained by geometrical operations on one site, which will be referred to as X1.
o X2 are obtained by inversion, X3 obtained by reflection in the x-y or x-z plane, and X4 by a combination of the two reflections.
o Y1, Y2, Y3, Y4, Z1, Z2, Z3, Z4 are obtained by cyclic permutations of the axes on the operations to generate the X family (e.g.switch x for z axis and repeat the above procedure).
Table S1 lists the relative coordinates of each of the dodecahedral sites in the cubic garnet unit cell.A .cif file of the unit cell with labeled sites is also available upon request.(3) Contributions to magnetic anisotropy from each site is found by summing cos  ⃗ ⃗ for each of its nearest magnetic neighbors, where  ⃗ is the unit vector between the two ions and  ⃗ is the magnetization of the neighboring ion.Summing for all RE sites we get the following equation for magnetic energy of a garnet with mixed A and B RE ions, relative to the unmixed A garnet: where C and C are constants and α are direction cosines.N are the number of B ions in the X, Y, and Z sites, and: where N is the number of formula units of (A1-xBx)3Fe5O12 in the crystal (x is the fraction of B ions in the crystal).
o (001) Sites fall into two groups, with 2/3 of them in the α group and 1/3 in the β group. α: It is interesting to note that for the (111) and (001) cases, the magnetic anisotropy is uniaxial with the unique axis perpendicular to the growth face.
(5) Concentration (x) dependence -In the (110) ordering scheme (as an example), A and B ions would have different "sticking coefficients" in each of the α, β, and γ sites.
Then, each site would contain some fraction of B ions: for the α sites and analogous expressions for the other sites.
At low concentrations, the energy equation can be simplified and approximated for weak site preference to depend on concentration quadratically, scaling as x(1-x).
The symmetry reduction of the unit cell for ordering on different growth surfaces is visualized in Fig. S1.As the symmetry of the growth surface is reduced from (111) → (110) → (112), the dodecahedral sites are split into 2, 3, and 4 degenerate site groups, as described above.The symmetry of the unit cells of resulting crystalline materials grown on these surfaces have lower symmetry than the unordered cubic cell (space group: Ia3 d).Considering a garnet film which includes all the symmetry variants of the atomic ordering, the space groups of these ordered cells are summarized in Table S2.
These site-ordered, symmetry-reduced structures are not conventional superlattices but they represent a 3D ordering based on the geometrical orientation of the 24 dodecahedral sites in the unit cell.The ordered structures are formed during the growth of these thin films by the preference of arriving RE ions to order into inequivalent dodecahedral sites due to steric factors.that the growth-induced anisotropy is robust up to high temperatures 4 , which indicates that the growth-induced order is retained to high temperatures.The order and the resulting anisotropy can be lost by a sufficiently high temperature anneal (e.g.>600˚C) to allow for RE diffusion within the structure.growth directions.

Supplementary Note 2: Strain calculation and additional high-resolution x-ray diffractograms
Reciprocal space mapping of the (642)+ reflection shows that the thickest film with the largest lattice mismatch (EuIG, 42 nm) is fully strained within the plane of the film, since the film and substrate peaks have the same qx value, as shown in Fig. S2(a).Therefore, it is assumed in all calculations that thinner films or those with less lattice mismatch will also be fully strained.The out-of-film plane spacing, related to 2θ of the (444) reflection (Fig. S2(b)), can be used to calculate the lattice parameter and strain of the film using expressions for a rhombohedral distortion.
The presence of Laue oscillations on the film peak as well as low spread in the rocking curve, as shown in Fig. S2(c), indicate high crystalline quality, low degree of mosaic spread, and planar top and bottom interfaces.Mixed garnet films of all compositions show similar rocking curve widths, indicating similar quality (Table S3).With the Globalfit software, the film peak is fitted assuming the film is fully relaxed, as shown in The GGG lattice parameter used in fitting is as = 1.2376 nm.We calculate film lattice parameter and strain geometrically (without literature values for stiffness) since we can assume that the films are fully strained to match the substrate in-plane lattice parameter.The lattice parameter can be used to determine the composition by the lever rule based on the strained lattice parameters of the end members, EuIG and TmIG.From the fitted thickness, the growth rates of the EuIG and TmIG are found to be 239 shots per nm and 614 shots per nm, respectively.Table S3 outlines the shot ratios, compositions, thicknesses, and strains of the REIG films, and the rocking curve widths.For the determination of Ms for films grown on GGG by magnetometry, it is well known that the paramagnetic signal from GGG must be subtracted, or the substrate must be thinned down to reduce this signal. 6Fig. S4 shows the VSM hysteresis loop for TmIG before and after linear background subtraction of each of the four tails (for a symmetric full loop).
The error associated with measuring the Ms of a given film can be estimated by the following contributions.Sample mounting error ( ± , units of emu/V) is estimated from the sample standard deviation of measured VSM calibration factors.Fitting error ( ± , units of V) is the standard error for the linear fit of the saturated branches of the hysteresis loops.The total error, ( ± , units of emu), is thus the combination of these two errors: This analysis is only applicable for easy axis loops (in the field range of 0.02 T).Beyond ~ 0.2 T, GGG has a non-linear background, so it is difficult to determine the hard axis saturation, which according to SMR measurements occurs for these garnets at most around 0.34 T. This is why only easy axis loops are reported.Table S4 reports the coercivities of the unpatterned films from the easy axis loops.Saturation field is extracted from SMR by numerically fitting the signal to a macrospin model, which determines the equilibrium magnetization direction for any applied field magnitude and direction 8 .Error is estimated to be the field step.The data and fit of a representative curve are shown in Fig. S5(c).
The magnetoelastic constants,   , were determined from SMR measurements of anisotropy on series of EuIG and TmIG films grown on (111) substrates with different lattice parameters, including GGG (as = 1.2376 nm), YSGG (as = 1.2426 nm), NGG (as = 1.2505 nm), GSGG (as = 1.2554 nm) 4,9 .Fig. S6(a)(b) shows the strain series of EuIG films.From the linear variation of anisotropy as a function of strain we find the value of   by linear regression of this equation: From this analysis,  for EuIG and TmIG are found to be (2.45 ± 0.8)×10 -7 and (-1.1 ± 0.1) ×10 -6 if  is taken to be that of YIG (766 GPa) 10 .The calculated magnetostriction coefficients are about a factor of five lower than the reported values for these materials.This could indicate a non-ideal stoichiometry or  could deviate from the bulk value for YIG.Now, assuming there are n such axes, all with the same polar angle  , with azimuthal angles evenly distributed, i.e.  = , the total energy is The second term cancels out, proven below using complex numbers.
But we can write the sum as a finite geometric series, and the numerator cancels out, since  = 1, thus the second sum in the energy equation is zero.
By the same logic, the geometric sequence term cancels out, thus we are left with  .The final energy expression for n axes becomes: The system has three important angles -1) the vertical axis which is the symmetry axis of the cone drawn out by the anisotropy directions (set to be the z-axis,  = 0), 2) the magnetization angle and 3) the anisotropy cone angle.Thus, both azimuthal angles are important and appear in the expression.
Considering the anisotropy cone angle fixed, we can write the energy: In the system, average energy is determined by ε/n, and is uniaxial in nature.
Description of the planes of order are shown in Fig. S7.

Supplementary Note 6: Scanning transmission electron microscopy (STEM)
It is essential to consider the role of the symmetry-allowed variants of the cation-ordering in order to analyse the structure and properties of the film.As described in Supplementary Note 5, for a (111) film an individual variant of the site order would yield a tilted anisotropy whereas the combination of the three variants produces PMA.Variants also affect the interpretation of XRD and STEM measurements.The (111) EuTmIG films follow the scheme of Fig. S8a,c.STEM along the [1 10] in-plane zone axis, Fig. S9, shows no elementally resolved site ordering because the sample consists of all three variants.The intensity in a STEM image results from averaging all the atoms in the column along the beam direction.Collecting information through all three variants averages out any compositional ordering through the thickness of the TEM lamella.Indeed, we can conclude that the correlation length of these variants is on the scale of a few unit cells or less.The lamella of the (111)-oriented film was ~10 nm thick, but we do not observe cation order, indicating that the size of the variants is smaller than 10 nm, and may be on the scale of one or a few unit cells (the unit cell size is ~1.2 nm).
This precludes visualizing the c site order in (111) garnet films unless we had a lamella that is thinner than the correlation distance of the site-order.For the (110)-oriented EuTmIG film, compositional analysis was carried out by EDS, analyzing the peaks of Eu, Tm, and Fe chosen to minimize overlap (Fig. S12).Then, non-linear principal component analysis was applied to reduce Poisson noise 11 .Lastly, the filtered image was slightly blurred with a Gaussian filter.Evidence of ordering can be seen even in the raw images (Fig. S12), but image processing makes the order immediately evident to the reader.Quantitative analysis of HAADF images of the film was also performed.The HAADF is sensitive to atomic number, Z, such that the averaged column intensities scale with approximately Z 2 and the number of atoms present in the column 12 .Tm peaks should be higher intensity than Eu, since ZTm = 69 and ZEu = 63.
First, a linear background subtraction was applied to reduce the intensity change due to the thickness gradient of the TEM lamella created during ion beam preparation.Then, atom columns were identified and fit using open-source code "Pycroscopy" 13 .Lastly, atom columns were masked with circles, and the intensity of each atom column was summed to minimize error in intensity measurement 14 .These atom columns were then binned into A, B, C, and D types, and averaged.Peak identification and an intensity histogram are shown in Fig. S13.
To examine the site occupancy of columns A,B,C and D, 1131 columns were identified.A and C columns both contain sites from the  group, but have different densities of atoms so different contain  and  group sites, with C columns higher intensity than the B columns.Recalling (Note 2) that there are twice as many  sites as +  sites, we expect for a site-ordered EuTmIG with x = 0.5 that Tm should occupy all the +  sites and ¼ of the  sites, and Eu should occupy ¾ of the  sites.The summed intensities were averaged for the two groups of atoms to give intensities that should be greater for the Tm-containing sites, as shown in Table S6.The analysis does not yield a statistically higher intensity for all the Tm-containing sites, which is inconsistent with the clear site-ordering observed in EDS.The discrepancy may be explained if there are cation vacancies (VRE) or Fe antisite defects (FeRE) in the RE columns: even a small amount of these point defects would give a larger effect on column intensity than the Z-contrast between the Tm and Eu.Indeed, when we look at the Fe EDS, we see a slight indication that excess Fe could be preferring β and γ sites over α sites.
We also compared the Eu and Tm distributions locally in the hexagonal rings of A and B columns.atoms and there would be a statistical fluctuation in the number of each atom in the columns, but the overall trend is clear and was also found in other sections of the image.
Concerning the role of defects, we did not observe any dislocations across the entire visible lamellae of (111) or (110)-oriented films.This is consistent with many other TEM investigations of epitaxial RE, Bi or Y garnets in our prior work 15,16 .Furthermore, the strain in the films does not relax even for thicknesses of 10s of nm, according to the RSM data, suggesting that dislocations do not form.Hence we do not believe dislocations play a major role in the site occupation or growth-induced anisotropy.Fe 2+ is another possible point defect arising from oxygen deficiency.We performed X-ray absorption spectroscopy on TbIG films in a prior work 16 which showed that amount of Fe 2+ was less than about 3%.Further, we would expect all the films to show similar amounts of Fe 2+ or other point defects such as vacancies because they Structural relaxation for the three garnets (mixed garnet and two end members) was performed, accounting for the energy, per-atom forces, and stress computed by spin-polarized, collinear calculations at each update of the atomic positions.Then, a final spin-polarized, non-collinear calculation was performed on this static relaxed structure to evaluate the energy with a specific magnetization orientation.Including energy changes due to magnetostriction would have required accounting during structural relaxation for the energy, per-atom forces, and stress that arise from orientation of the magnetization.Such a method for structural relaxation was deemed prohibitively expensive since it would entail minimizing the energy and ensuring per-atom forces approach zero with spin-polarized, non-collinear electronic structure calculations performed at each update of the atomic positions rather than merely at the final step at the end of the structural relaxation.As a result, the calculation method presented does not include this magnetostriction and elastic energies.However, the mixed garnet has a smaller magnetostriction coefficient than the end members, and magnetostriction would therefore be unlikely to account for the larger anisotropy energy observed for the mixed garnet by DFT.Fig. S15 shows that for a site-ordered unit cells with fixed RE composition (Eu1.5Tm1.5Fe5O12),as the degree of ordering of the RE on distinguishable sites increases, the intensity of the (110) order peak also increases with a quadratic dependence.The peak intensity depends on the structure factor which varies with the difference in atomic number of the ordered RE cations.For the completely mixed case, there is no difference in the average atomic number on each site, so no peak is seen.Lastly, to explain the observation of weak (110) peaks in the asymmetric scans of endmember EuIG and TmIG films grown on (111) GGG (Fig. 3), we consider the effects of point defects.
Generally, removing atoms from a unit cell reduces the symmetry of the cell, resulting in the appearance of structurally forbidden peaks.Fig. S16 shows the consequence of removing a single oxygen atom from the unit cell of GGG, for example.Even the removal of one atom in 160 produces sufficient change in the structure factors to produce peaks which were absent in the perfect crystal.
We do not expect many defects in the Czochralski-grown GGG substrates, and the (110) peak of the GGG in Fig. 3 is absent.However, YIG or REIG films growth by pulsed laser deposition often exhibit non-ideal cation stoichiometry or oxygen deficiency. 16The presence of a weak (110) peak for the endmembers EuIG and TmIG is readily explained by such point defects.Point defects are also expected to contribute to the (110) peaks in the EuTmIG film, but the much greater intensity of the peak is attributed to the Eu/Tm site ordering.For the GGG substrate and the films, (220) and (440) reflections are present as expected since these even reflections are not systematically absent.The (110) peak (at 2θ = ~10°) is normally forbidden for the Ia3 d cubic garnet crystal structure, but it is present for all samples.This is due to the well-known phenomenon of umweganregung, which causes the appearance of normally forbidden peaks in symmetric scans due to reasons other than symmetry reduction in the crystal (Fig. S17b). 19The presence of the umweganregung peak prohibits us from using the (110) reflection to verify cation order in these films, since it is also present in the substrate.[It is important to note, however, that the (110) peak can still be used to diagnose order in the (111) films as shown in Fig. 2 because it is collected in a skew geometry, which reduces the umweganregung peak.] In contrast, the higher order reflection, (330), is much less prone to umweganregung.This peak is therefore a good diagnostic of cation order.We note that (330) peaks exist for the EuIG and EuTmIG films on GGG, but not the uncoated GGG substrate (Fig. S17c).Moreover, the dspacings of these (330) peaks are consistent with the higher lattice parameter of the EuIG compared to the EuTmIG, unlike the (110) peaks, giving us confidence that the (330) peaks are not umweganregung peaks from the substrate.The (330) peak of the EuTmIG is significantly higher than that of the EuIG.We believe that the EuIG shows some (330) intensity due to point defects (Supplementary Note 8), but the higher intensity of the EuTmIG peak is a result of RE site ordering.Thus, we can verify cation site order in the (110)-oriented EuTmIG from both XRD and STEM.

Note 1 : 2 Note 2 : 8 Note 3 : 4 : 5 : 6 :Note 7 : 8 . 9 .
Callen's theory of growth-induced anisotropy (GIA) Strain calculation and additional high-resolution x-ray diffractograms Vibrating sample magnetometry background subtraction and error propagation Note Spin Hall magnetoresistance and anisotropy calculations Note Derivation of uniaxial anisotropy from three tilted anisotropies Note Scanning transmission electron microscopy (STEM) Density Functional Theory (DFT) Calculations Note Simulations of X-ray Diffraction (XRD) from ordered garnets Note Relationship between XRD and STEM data for verification of cation order References Supplementary Note 1: Callen's theory of growth-induced anisotropy (GIA)

( 4 )
For different growth faces, c sites fall into symmetrically inequivalent categories according to how the site symmetries are reduced at the surface.The anisotropy energy can be simplified for each case in which a B ion is substituted in an A site.The c sites are reduced into groups α, β, γ, δ of degenerate sites according to the symmetry of the neighbors around each site at the growth surface.The sets of equivalent sites and the simplified expressions for energy for several growth surfaces are: o (110) Sites fall into three groups, with 2/3 of them in the α group and 1/6 in each of the β and γ groups.α: NX1 = NX3 = NY2 = NY4, NX2 = NX4 = NY1 = NY3 (5)

Fig. S1 .
Fig. S1.Order schemes and inequivalent dodecahedral site groups for the [111], [110], and [112] Fringes also allow us to fit thickness of the film, analogously to Kiessig fringes in thin film reflectivity.Fig.S2(d)shows that only garnet peaks corresponding to the out of plane lattice direction are present over a wide range of angles, confirming the phase purity of the films.

Fig. S3(
Fig. S3(a), although reciprocal space mapping shows that films are fully strained.This fit gives us the relaxed lattice parameter, and hence d , the lattice spacing in the out of plane direction of the film (Fig. S3(b)):

Fig. S3 .
Fig. S3.Geometric representation of a strained (111) unit cell on a substrate.(a) Schematic of fully relaxed and fully strained films with the same out of plane lattice spacing.(b) Geometric relation between out of plane spacing and relaxed lattice parameter.(c) Geometric relation for in plane spacing, d.(d) Rhombohedral relation for strained film.

Fig. S7 .
Fig. S7.Views of the individual tilted anisotropy axes for the variants of the ordered (111) REIG system (top and side views).
Fig. S8 illustrates variants for the case (111) and (110) films.In (a) an ordered unit cell of one variant of the (111) case is shown with a reference plane in red; red and blue balls represent the inequivalent c-sites.In Fig. S8(c), multiple unit cells are shown for the (111) ordering, each corresponding to one of the three variants.Because of the superposition of the three variants, there is no global order in the c-sites that can be measured by EDS unless the spatial extent of the variants is large compared to the sample.In contrast, for the case of (110) order, Fig. S8(b) shows one variant in which red, blue and purple balls represent the three sets of inequivalent c-sites .Fig. S8(d) shows a film consisting of the two variants.If we image along one specific zone axis it is possible to identify the order by EDS: we can distinguish red from blue/purple sites, since long range ordering of the red site motif prevails, but blue sites cannot be distinguished from purple.

Fig. S10 .
Fig. S10.[11 1] zone axis of (110) grown EuTmIG on GGG.(a) View of interface with 10 nm scale bar.(b) View of interface with 5 nm scale bar.(c) View of interface with 2 nm scale bar.In each figure, the interface appears darker.

Fig. S12 .
Fig. S12.Unprocessed EDS of Eu, Tm, and Fe in the ordered EuTmIG film and total EDS spectra.

Fig. S13 .
Fig. S13.Atom column intensity comparison.(a) selected atoms from a single image (b) histogram of atom column intensities for A, B, C, and D sites.(c) EDS of the (110) EuTmIG, showing site preference for Eu on A sites and Tm on B sites.(d) Extracted intensity line scan for Eu and Tm along the path of the hexagon in (c).
Fig.S13(c,d) shows the analysis of one region of the sample with a hexagon of dodecahedral sites.Columns 1, 2, 4 and 5 contain sites from the  group and columns 3 and 6 contain sites from the  group.The line plot clearly shows that columns 3 and 6 contain the least Eu and most Tm, and the other columns have more Eu and less Tm.The peaks are not exactly the same for columns 3 and 6 (or for columns 1, 2, 4 and 5) because the column contains ~20

Fig. S16 .
Fig. S16.Simulated powder diffraction of gadolinium gallium garnet in the perfect state (upper panel) and with one oxygen atom removed from the unit cell, forming a vacancy (lower panel).

Table S1 .
Fractional coordinates of labelled dodecahedral sites in the garnet unit cell

Table S2 .
Crystallographic space groups of cation-ordered garnets.

Table S3 .
Summary of structural properties for REIG films

Table S6 .
Averaged atom column intensity and relevant statistics.